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A real-time path integral Monte Carlo approach is developed to study the dynamics in a many- 
body quantum system until reaching a nonequilibrium stationary state. The approach is based on 
augmenting an exact reduced equation for the evolution of the system in the interaction picture 
which is amenable to an efficient path integral (worldline) Monte Carlo approach. Results obtained 
for a model of inelastic tunneling spectroscopy reveal the applicability of the approach to a wide 
range of physically important regimes, including high (classical) and low (quantum) temperatures, 
and weak (perturbative) and strong electron-phonon couplings. 



In recent years there has been considerable interest in 
the study of quantum mechanical systems at the nanome- 
ter scale that are driven out of equilibrium. Experimental 
breakthroughs on transport in molecular junctions have 
uncovered fascinating behavior in molecular systems far 
from equilibrium [1, 2J. Much attention has been paid 
to the study of the transport through strongly correlated 
systems in which electron correlations are dominant and 
lead to interesting physics such as the nonequilibrium 
Kondo effect [3 HI [51 |B] and Coulomb blockade |7j in 
quantum dots/molecules, tunneling in a Luttinger liq- 
uid [S] [3] , or inelastic effects induced by electron-phonon 
interactions |10j as probed by inelastic electron tunneling 
spectroscopy [TT]. 

Exact theoretical treatment of such many-body sys- 
tems is rather sparse and includes only a small class of 
simplified model problems [HJ [T31 [S]. For a general 
solution mean field equations based on the many-body 
nonequilibrium Green's function (NEGF) approach can 
be formulated [151 and solved numerically. However, the 
mean-field NEGF approach is quite limited since in many 
cases it is based on a perturbative scheme and the inclu- 
sion of higher order corrections is not always clear or sys- 
tematic. Traditional determinant Monte Carlo (MC) ap- 
proaches based on standard discretized Pis [111 IlZj have 
been used but are numerically rather expensive and seem 
to be restriced to imaginary-time calculations. Thus, the 
development of a general approach suitable for the treat- 
ment of nonequilibrium many-body quantum systems re- 
mains a grand challenge. 

In this letter we present a novel approach aimed at 
obtaining exact numerical results for various dynam- 
ical quantities such as the current I{t), conductance 
g(t), dot population P(t), etc. in a many-body quan- 
tum system that is driven out of equilibrium. We fo- 
cus on a well studied model of inelastic tunneling spec- 
troscopy mi nil [201 [21] I where a quantum dot is coupled 
to two fermionic reservoirs representing the left and right 
leads at chemical potentials of /i^ and /i/j, respectively, 
and to a bosonic bath representing the phonon environ- 
ment. Motivated by the success of real time PIMC tech- 
niques developed for molecular chains [22j . we adopt a 



similar procedure and formulate an exact real time path 
integral (PI) representation for the dynamical quantity 
of interest. To reduce the computational complexity we 
integrate out the fermionic leads and the bosonic envi- 
ronment and obtain expressions for their corresponding 
influence functionals [23] [23] . We develop an adequate 
Monte Carlo procedure where we propagate the density 
matrix from an initial factorized arbitrary condition to- 
wards a steady-state. 

A nonequilibrium quantum dot in a phonon environ- 
ment can be described by the Hamiltonian 

= ^ heua\ak+ ^ {tt,a\d + \i.c^ + hen (i^ d 
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The time-dependent current from the left lead onto the 
dot can be written as 

Idt) = -ej^{NL{t)) - -|lm^t,.(4(t)d(t)) 

keL 

= -yIm^tfce"'=*tr{M/o;7j(t)4dff,(t)C//(i)} ,(2) 
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where Wa = W^l*^ x '^D^Ph ^^e initial factorizing prepa- 
ration, Ui{t) ~ giHat/h^-iHt/h jg ^Yic interaction picture 
time evolution operator with Hq = i?LR + -ffD,Ph = 
ffLR+i^D+i^^''ph+■^Ph, and dH„{t) = e'Hot/n^^-iHot/n^ 
We present the approach for a quantum dot which is 
initially empty; the expressions for the case of an initially 
occupied dot (or a mixed preparation) can be obtained 
straightforwardly, as well as those for the right current, 
the conductivity, the dot's population etc. After express- 
ing the time evolution operators in Eq. ([2| by virtue of 
Dyson series, I Lit) can be written as an infinite sum over 
time-ordered integrals. 
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Here, dn^^^it) = e^H^-^^*/^de-'"°-^'-*''^ = dH„{t). The 
trace over the leads degrees of freedom can now be per- 
formed exactly (see, e.g., Ref. US]), yielding 



C{ti,...,t2N) det(Af) 
where M is a matrix with elements Mij — 
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t2i) + ^i{t2,-l-t2^){l-St,^_,,t){^-^t,,,t)] for J < Z < 
or [S>(^2,-l-^2^) + S>(^2J-l-^2^)(l-<5^,,„l.t)(l-<5t2..t)] 

for j > i. Here, E^^(t) (S^j^(t)) denotes the leads' 
lesser (greater) self energy in the time domain. Q in 
Eq. ([3| represents a 2(n-|-n'-|-l) point correlation function 
of the dot-phonon subsystem along the Kadanoff-Baym 
contour s : — *■ i — *■ 0. It is most conveniently evalu- 
ated in the framework of Feynman Pis, which allows to 
integrate out the phonon degrees of freedom exactly |23j , 
yielding 



g = exp|ie£i J ds [a{s) - a'{s)] | J^[a, a'] , 



(5) 



where denotes the Feynman- Vernon influence func- 
tional |23j summarizing the influence of the phonons on 
the dot, and a, a' e 0, 1 denote the corresponding for- 
ward and backward branches of the dot path, referring 
to the propagations s : — > t and s' : t 0, respectively. 
Since -ffD,Ph lacks any terms which could flip the state 
of the dot from empty (cr = 0) to occupied (cr = 1) or 
vice versa, the contour-ordered sequence of time points 
si, . . . , S2n+i,t, S2n'i • ■ • J '^1 Uniquely deflnes the paths a 
and cr': Every refers to one occurrence of the dot- 
state altering interaction part H^^]^^ in the Dyson series. 
Therefore, the {sj, s'j,} define the kink times of the dot 
path pB] , between which (t(s) remains constant. In this 
spirit, Eq. (|3| closely resembles the instanton expansion 
of the partition function in the spin-boson model [27^ . 

While the expressions Q and ^ for £ and Q allow 
for a rather compact expression of iL^t), they introduce 
retardation effects which are arbitrarily long ranged in 
time, making analytical progress rather cumbersome. To 
allow for a numerical evaluation of Eq. ^ , the integrals 



in Eq. ^ are approximated dividing the time axis into 
q discrete steps, 

/ dsN / dsN-1--- / dsif{si,. . . ,sn) 
Jo Jo Jo 

where r = t/q. While this introduces a systematic er- 
ror of the order of 0(t) and establishes an upper bound 
to the sum over the kink numbers, the corresponding 
errors can be made arbitrarily small by adjusting q to 
correspondingly large numbers. In addition, systematic 
improvement can be made by using a more accurate inte- 
gration scheme. Performing the sums over all kink num- 
bers and the corresponding (discretized) kink times is 
now equivalent to summing over all possible (discretized) 
dot paths {cjj — a{jT),a'-, = cr'(j'r)}, such that Eq. ([s" 
can be written as a discretized PL 
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where rikink denotes the number of kinks of the path 
{aj, a-'ji}- Eq. (ITj) can now readily be evaluated by means 
of PI (or worldlme) MC [22l[26]. 

We now turn to discuss the application of the pro- 
posed approach to the model system described by the 
Hamiltonian ([ij. In the present approach the effect of 

the leads is fully determined by the self energies S^^^(i) 
(cf. Eq. Q), which are defined in terms of r(e) — 
+ rfl(e) in Fourier space [TS^. r(e) is taken to be 
energy independent (wide band limit) with a soft cutoff 

at e = ±6^: TLfR{e) = [i^,.(.-^c')n(+e-^(^+^e)] ■ In all re- 
sults presented below, F^ = F^, ly = 5T = 5(Fl + Fj^), 
and Cc — lOF or 20F to converge the results. Simi- 
larly, the phonon influence functional is completely 
specifled by the phonon spectral density [53], J{uj) — 
W M^S{uj — LUa). For a single phonon coupled to a 
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FIG. 1: Plots of the time dependent current /^(t) (lower left 
panel), (lower right panel), average current I{t) (upper 

left panel), and the dot population P{t) (upper right panel) 
for Mo = 4 and fesT = ^, in units of T = + Tr. 



secondary bath via a coupling constant 7, J{lo) becomes 
a Lorentzian: J(llj) — 7^ ^r? 

^ ' [(c^/a;o)2-l]2 + [fi,2^,^oc^/(2M2)]' >— ' 

We note in passing that the proposed simulation scheme 
does neither depend on the particular form of E^^^(t) 
nor of J{uj). 

In Fig. [1] we plot the time dependent current and the 
time dependent dot population for different values of 
the model parameters for T < T (quantum regime) and 
Mq > r (strong coupling). Lower panels show the left 
{I hit)) and right (Ij^(t)) currents and upper panels show 
the average current I{t) = ^[lL{t) + Init)] and the dot 
population P{t). The left/right currents are character- 
ized by damped coherent oscillation with a long time ex- 
ponential decay to a steady state value. The present MC 
scheme provides converged results for the time-dependent 
current despite the notorious dynamical and fermionic 
sign problem. This can be contributed to electronic de- 
phasing induced by the leads and the bosonic bath, as 
well as to the rather small size of the determinants in 
Eq. ^ (i.e. half the number of kinks on the combined 
forward-backward dot path) which allows for very fast 
MC sampling. While the steady state current can be ex- 
tracted from an exponential fit to /L/ij(i), the average 
current I{t) typically decays much faster to steady state, 
such that in most cases the steady state value could be 
obtained as the corresponding plateau value. 

The fact that I{t) approaches faster to steady state is 
consistent with other flux-based methods, and can 
be rationalized by looking at dynamical fluctuations un- 
der equilibrium conditions. This is depicted for the case 
where fJ-L = fJ'R and the system decays from an initial 
factorized density to equilibrium. The left/right currents 
show pronounced coherent oscillations with flnite values 
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FIG. 2: Plots of the total current 7 as a function of the bias 
voltage for fieo — 0, A*i,fl ~ ^1 ^'iid ksT = |, in units 

of r = Pi + Tr. 



even at i > bh/T while their average value vanishes for 
all times as appropriate for equilibrium conditions. 

To analyze the limitations of the approach we have 
conducted a systematic study for the case where Mq — 
(not shown), which numerically is the most difficult case 
since decoherence effects arising from electron-phonon 
coupling are neglected. Comparison of the numerical re- 
sults and the corresponding analytical solution reveals 
that the present PIMC scheme provides converged re- 
sults for a wide range of gate and bias voltages, and for 
a wide range of temperatures spanning the classical limit 
down to the quantum regime. The method is accurate as 
long as I{t) decays exponentially to steady state within 
the time window of i < tmax- As can be seen in Fig. [T] 
the decay of I{t) is slower for decreasing values of wq 
and fJ-L — f^R- The approach is still limited for small val- 
ues of these parameters, as well as other regimes where 
the steady state timescale and/or the decay time of co- 
herent oscillations are similarily stretched (e.g. very low 
temperatures); however, improvement of the MC scheme 
should, in principle, overcome these limitations. 

In Fig. [2] we plot the steady state current / as a func- 
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tion of the bias — Mi?, = eV for different electron- 
phonon couplings Mq (lower panel) , different phonon fre- 
quencies ujq (middle panel), and for different couplings 7 
to a secondary bath (upper panel) for ksT < T (quan- 
tum regime) . We compare the results to an approximate 
method based on a generalization of the single particle 
approximation |30j to include the leading order term of 
the Fermi sea [TB]- 

When the coupling between the primary oscillator and 
the secondary phonon bath is small we observe steps 
in the voltage dependent current at integer values of 
eV = 2H0J0 (middle panel) [TD]. As the coupling to the 
secondary phonon bath increases or for an Ohmic spec- 
tral density the steps diminish and eventually disappear 
(upper panel), signifying the wide range of phonon fre- 
quencies contained in the spectral density. The results 
for the steady state current shown in the lower panel of 
Fig. |2] span a wide range of electron- phonon couplings 
from the Landauer inelastic case, through the perturba- 
tive regime to the strong coupling limit. As Mq increases 
the value of the current decreases from the Landauer in- 
elastic single-channel value to lower values. Our approach 
clearly captures elastic effects at all values of Mq > as 
depicted by the lower values of the current and by the 
steps at twice the frequency of the primary phonon mode. 

Comparing the exact numerical real-time PIMC results 
to the approximate expansion method we find that agree- 
ment is quantitative for small values of the voltage cor- 
responding to the first step of I{V). For larger values we 
observe significant deviations between the exact numer- 
ical results and the approximate method for Mq/F > 4. 
These deviations signify the importance of high order ef- 
fects of the Fermi sea to the transport process, which are 
naturally included in the real-time PIMC approach. 

In conclusion, we have developed a novel real-time 
PIMC approach to study the dynamics in open quantum 
systems that are driven out of equilibrium. The approach 
is based on expressing the time evolution by virtue of 
Dyson series, before reducing the dynamics of the entire 



system by integrating out the fermionic/bosonic baths 
and introducing exact influence functionals. The remain- 
ing infinite sum over contour-ordered time integrals is 
then evaluated by PI (worldline) MC. We have applied 
the approach to study the time-dependent current in a 
well-studied model of inelastic tunneling spectroscopy, 
where a quantum dot is coupled to two fermionic leads 
and to a bosonic phonon bath. Numerical results indi- 
cate that the approach is robust and can be used for 
a wide range of model parameters spanning the classi- 
cal to quantum limits, a range of experimentally acces- 
sible chemical bias, different phonon frequencies, weak 
to strong electron-phonon couplings, and a wide range 
of couplings between the primary phonon mode and a 
secondary phonon bath. Furthermore, the approach pro- 
vides real-time information for various quantities, includ- 
ing the current, conductance, electronic population, spec- 
tral function, and more. 

We believe that our approach is capable of resolving 
several shortcomings found in currently used approaches. 
First, it is not based on any perturbative treatment and 
in principle can provide exact numerical results. Second, 
it yields a compact expression for the dynamics which 
can be evaluated by the proposed PI worldline method 
in a numerically very efficient way. Furthermore, the 
real-time propagation scheme allows to study transient 
phenomena and timescales and also to include time de- 
pendent fields. In addition, since it is based on a MC 
procedure, an enhancement of the accuracy can be ob- 
tained by improving the sampling scheme, as well as by 
including already existing schemes to sooth the dynami- 
cal and fermionic sign problem. Finally, it can be applied 
to a general many body problem, as long as a stable PI 
worldline approach can be derived. 
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